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ABSTRACT 

We study the observable signatures of self-gravitating MHD turbulence by applying the probability 
density functions (PDFs) and the spatial density power spectrum to synthetic column density maps. 

We find that there exists three characterizable stages of the evolution of the collapsing cloud which we 
term “early,” “intermediate,” and “advanced.” At early times, i.e. t < 0.15t//, the column density has 
a power spectral slope similar to nongravitating supersonic turbulence and a lognormal distribution. 

At an intermediate stage, i.e. 0.15t// < t < 0.35t//, there exists signatures of the prestellar cores in 
the shallower PDF and power spectrum power law slopes. The column density PDF power law tails at 
these times have line of sight averaged slopes ranging from -2.5 to -1.5 with shallower values belonging 
to simulations with lower magnetic field strength. The density power spectrum slope becomes shallow 
and can be characterized by P{k) = where Ai describes the amplitude, describes the 

classical power law behavior and the scale kc characterizes the turn over from turbulence dominated to 
self-gravity dominated. At advanced stages of collapse, i.e. ~ t > 0.35t//, the power spectral slope is 
positive valued, and a dramatic increase is observed in the PDF moments and the Tsallis incremental 
PDF parameters, which gives rise to deviations between PDF-sonic Mach number relations. Finally, 
we show that the imprint of gravity on the density power spectrum can be replicated in non-gravitating 
turbulence by introducing a delta-function with amplitude equivalent to the maximum valued point 
in a given self-gravitating map. We find that the turbulence power spectrum restored through spatial 
filtering of the high density material. 

Subject headings: methods: numerical — AMR, MHD 


1. INTRODUCTION 

Molecular clouds are highly turbulent, magnetized 
and are the sites of all known star formation 
(|Elmegreen &: Scal3l2004fl . The details of the collapse of 
molecular clouds determine the key properties of the star 
formation rate (SFR) and stellar initial mass d istribution 
(IMF, see e.g. iHennebelle fc Falgarond 1201^ 1 Thus, the 
development of a detailed understanding of the dynamics 
of molecular clouds is an essential step toward a complete 
picture of star formation, including predicting the initial 
mass function. 

The turbulent nature of molecular clouds is ev¬ 
ident from a variety of observations including 
non-thermal broadening of molecular emission and 
absorption li nes such as carbon monoxide (see 


Snitzed 119781 IStutzki & Guested 119901: iHever & BruntI 

20041! and fractal and hierarchical structures (see 

Elmegreen & Ehnegreed 119831 IVazauez-Semadenil 11994 

Burkhart et al.l 120131! A 

number of new techniques, 
the turbulence velocity 
OOOl for a review) and the 
Alfven Mach number (see 

including those studvinj 
soectrum fsee iLazariad 12 
sonic Mach number and 

Kowal & Lazariad 120071 

Burkhart et al.l 120091 120101 

20121 lEsauivel & Lazariad 

20101: ITofflemire et aLII201 


have been applied to Giant Molecular Clouds (GMCs) 
and also have shown that turbulence there is supersonic. 
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In light of this, turbulence seems to play a duel role 
in the CMC environment of providing support on the 
large scales (i.e. scales of t he cloud) while compressiir g 
small scales via shocks (IMac Low &: KlessenI 120041! . 
Althoug h it is clear that molec u lar clouds are tur- 


bulent (lElmeffreen & Scald 

20041: IMac Low & Klessen 

2004D. mametized (iCnitcher 

20121! and self-gravitating 

(ILarsonl I1981I: ISolomon et alJ 11987ft. the relative 
importance of these components is still under de- 

bate (e.g. IPadoan & Nordlundl 119991: ILi et al.l 

2009; 

VazQuez-Semadeni et al.l 120081: iKritsuk et al.l 

2013; 


Li et al.l 120141! . As a mature theory of magnetized, su¬ 


personic, self-gravitating turbulence has yet to emerge, 
interpretations of observations are difficult. In light of 
this, numerical models have proven to be an important 
tool in the understanding of observations. 

Several recent numerical studies have probed the 
observational signatures of turbulent clouds, both with 
and without self-gravity and with a range of Mach 
numbers and magnetic field strengths. Simulations, 
observations and theoretic works have all pointed out 
the fact that compressible turbulence is important 
for creatin g filaments and re gions o f high density 
contrast (IKow al fc La zarianI l2007t iBurkhart et al.l 
l2009t iFederrath et al.l 120101! . Shocks can broaden the 
density/column density probability density function 
(pdf) and shallow the power spectral slope. The 
column density PDF and power spectrum (or delta 
variance, which is isomorphic to the power spectrum) 
are the chief targets for this study, due to their im¬ 
portance to tu r bulence based star format io n theories 
(iPadoanI 119951 : iKrumholz fc McK^ 120051 : iFederratB 
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l2013f) . Both the density PDF and power spectrum 
have been shown to be sensitive to Ma ch number, mag¬ 
netic field strength, and self gravity (iBeresnvak et alJ 
2005t_^owal&jjazari^|2007|jKainulainei^£^aIJ 2009; 


Burkhart,etM ■I l20ir)t iBurkhart fc LazarianI l2012t 

Collins et alJl2012t iFederrath fc KlessenI 120131 1 however 


less attention has been paid to higher order statistics of 
column density. 

In this work we explore the variation of the properties 
of the PDFs and the power spectrum of column density as 
simulated clouds collapse under their own self-gravity as 
a function of time and magnetic field strength. The goal 
of this paper is to provide observers with useful meth¬ 
ods to apply to dust extinction maps and integrated in¬ 
tensity maps, and to explore if the PDFs/spectrum of 
column density behave si milarly to 3D densit y studied 
in previous works such as iCollins et al.l ()2012[ 1 . In par¬ 
ticular we are interested to explore if the collapse evo¬ 
lution of the cloud can be observationally discerned and 
pay particular attention to the time dependency of the 
turbulence statistics. Further, we explore two novel fits 
to the column density power spectrum, and use one of 
these techniques to disentangle the signatures of self¬ 
gravity from the turbulence and to shed light on the ori- 
gin of the increase in the d ensity power spectrum slope 
(|Federrath fc KlessenI I2013D . The ultimate aim of this 
paper is to provide set of tools to quantify the physical 
state of observed molecular clouds. 

The layout of the paper is as follows. In Section [2] we 
describe the simulations studied and the codes used to 
generate them. In Section[3]we introduce the the column 
density PDF and apply it to our simulations. We inves¬ 
tigate the evolution of the PDF power law tail (Section 
1^ . the evolution of the variance (Section 13. 3|) , higher 
order moments ^Section 13. 4|) and finally the Tsallis dis¬ 
tribution (Section [331). Section [4] we discuss the evo¬ 
lution of the slope of the power spectrum, a new fit to 
the evolving spectrum fSection l4.ip : a technique to re¬ 
produce the self-gravitating spectrum from the turbulent 
spectrum ISection 14.211 : and a technique to recover the 
turbulent spectrum from the self-gravitating spectrum 
. We discuss the implications of our results 


((Section 14.3 

in Section |S] followed by our conclusions in Section [51 


Alfvenic (i?ext = 1-0) and super-Alfvenic (Bext = 0.1) 
turbulence. 

The three self-gravitating simulations were performed 
using the c onstrained transport MHD option in Enzo 
(MHDCT) ((Collins et al.ll20inl: iBT^n et ahllMl . This 
co de uses the sec ond order hyperboli c solver described 
bviLi et al.l (|2008ll . the CT method of iGardiner &: Stolii 
(1200,511. and th e divergence free interpolation method of 
iBalsard (1200111. The s e simu lations were described in de- 
tail in ICollins et al.l (|2012l l. The self-gravitating Enzo 
models are supersonic and super-Alfvenic, with mass 
chosen such that kinetic and gravitational energies are 
equal. They have three different values of initial mag¬ 
netic field which sets the Alfven Mach number. Through¬ 
out the text and in Table 1, we denote these different 
cases as low, mid, and high to stand for “low-valued,” 
“middle-valued,” and “high-valued” magnetic field runs. 

In all simulations, the ideal MHD equations are solved 
in a periodic box, using an isothermal equation of state 
(p = c^p) and a variety of sonic and Alfvenic Mach num¬ 
bers (Ms = (|v|)/(cs) and Ma = i/p(|v|)/(|B|), respec¬ 
tively). Here, p is density, v is velocity, B is magnetic 
field, p is the gas pressure, and Cg is the isothermal speed 
of sound. In all simulations, periodic cubes with ini¬ 
tially uniform density and magnetic fields are driven with 
solenoidal forcing, using the ideal MHD equations. The 
self-gravitating simulations were driven until a steady 
state was reached, at which point gravity was turned on. 
Driving continued during the collapse phase. A summary 
of the simulations can be found in Table |T] Four repre¬ 
sentative snapshots of the Enzo simulations can be seen 
in Figured) 

For the Enzo simulations we select the Mach num¬ 
ber, Ms-, virial parameter, Ovir, and mean thermal-to- 
magnetic pressure ratio, /3o as 


Ms 

Ovir 


/3o 


= 9 


bvz 


3GpoAq 

Sttc^Po 


= 1 


= 0 . 2 , 2 , 20 , 


( 1 ) 

( 2 ) 

(3) 


2. METHOD 

Two suites of simulations were used for this paper: one 
suite of non-gravitating simulations that used a fixed res¬ 
olution of 512^ (which we refer to as the Godunov simu¬ 
lations), and one suite of self-gravitating simulations that 
used a 512^ root grid and four levels of adaptive mesh 
refinement (AMR, which we refer to as the Enzo simula¬ 
tions). All simulations solved the ideal MHD equations 
with large-scale solenoidal forcing. 


The Godunov simulations used here are ten non¬ 
gravitating simulations with sonic Mach numbers ~ 
0.5 — 22 and have b e en u se d in many p r evious 


works (ICho & Lazarian 

120031: IBurkhart et al.l I2009L 

120101 iKowal & LazarianI 

20071: IKowal et all 120091 EfTill. 


These simulations us ed the algorithm described in 
(|Gho fc Lazarianl2002t) , which uses a combination of spa¬ 
tially third-order essentially nonoscillatory (ENO) meth¬ 
ods to ensure both accuracy and stability, and a third- 
order Runga-Kutta time integration. The Godunov mod¬ 
els are divided into two groups corresponding to sub- 


where Wrms is the rms velocity fluctuation, po is the mean 
density, Lq is the size of the box, and Bq is the mean 
magnetic field. These can be scaled to physical clouds 


as 

tg = l.ln“_Y^Myr (4) 

Lo = 4.6cs,2?T-H,V^PC (5) 

Vrms = 1.8cs,2km s“^ (6) 

M = 5900cs.2nH,V^MQ (7) 

Ho = (13,4.4,1.3)cs.2nH(3pG, (8) 


where Cs ,2 = 0.2km s“^ and nn.a = n///(1000cm“^) are 
the sound speed and hydrogen number density, respec¬ 
tively, and we have used a mean molecular weight of 2.3 
amu per particle. 

In what follows we use column density maps from the 
Enzo and Godunov simulations to investigate the util¬ 
ity of often used statistical descriptors of turbulence and 
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Table 1 

A summary of the simulations presented here. 


Model 

Ms 

Ma 

Description 

1 

0.5 

0.7 

subsonic & sub-Altvenic, Godunov 

2 

4.2 

0.7 

supersonic sub-Alfvenic, Godunov 

3 

8.0 

0.7 

supersonic sub-Alfvenic, Godunov 

4 

10 

0.7 

supersonic &; sub-Alfvenic, Godunov 

5 

20 

0.7 

supersonic sub-Alfvenic, Godunov 

6 

0.5 

2.0 

subsonic & super-Alfvenic, Godunov 

7 

4.2 

2.0 

supersonic Sz super-Alfvenic, Godunov 

8 

8.0 

2.0 

supersonic & super-Alfvenic, Godunov 

9 

10 

2.0 

supersonic Sz super-Alfvenic , Godunov 

10 

20 

2.0 

supersonic Sz super-Alfvenic, Godunov 

high 

8.5 

6.5 

supersonic & super-Alfvenic, Enzo 

mid 

8.2 

12. 

supersonic & super-Alfvenic, Enzo 

low 

9.1 

30 

supersonic & super-Alfvenic, Enzo 



Figure 1. Column density plots for the high magnetic field Enzo simulation. Time is shown in the bottom left. 
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gravitational collapse, including the PDFs, Tsallis dis¬ 
tribution and Fourier power spectrum, for observations 
of column density, e.g. dust extinction maps. We ac¬ 
count for fluctuations along different sight-lines relative 
to the orientation of the mean field in error bars, taken 
as the standard deviation between a measure along dif¬ 
ferent sight-lines. As observers often do not have in¬ 
formation regarding the mean field we treat the line of 
sight relative to the mean magnetic field as an unknown 
parameter. Regarding the ENZO simulations, they are 
all super-Alfvenic and thus any anisotr opy introduced 
by the mean field will be very weak (see iBurkhart et al.l 
(|2014f) l. The Godonov simulations do have sub-Alfvenic 
snapshots, however the effec t of anisotropy in these is 
studied in other wor ks (see lEsauivel fc LazarianI 120111 : 
IBurkhart et al.ir2ni4l ll. 

3. PDFS 


The Probability Density Eunction of gas and dust in 
star forming regions provides important signatures of 
both MHD turbulence and gravitational collapse. In 
the case of low density MHD turbulence, where self- 
gravitational influences are negligible, the PDE exhibits 
a lognormal form for both the density and column den¬ 
sity distributions, i.e. 

= 7i=f 


where ^ = — a^ /2 is the mean of In p and a is the stan¬ 
dard deviation (jBlaisdell et al.lllQ^ IVazauez-Semadeiiil 

[HU). 

Lognormal column density PDFs have been observed 
in various phases of the ISM including in molecular 
(IVazauez -SemadenillI994 lBruntll20 I0l: IKainu lainen e t ahl 


2014^ I2012t IMolina et al.l 120121 : 

Kainulainen &: Tanll2013h and in the di f fuse warm neu- 


tral and ionized ISM (|llill et al.1120081: IBurkhart et al.1 
(Moll . At the onset of gravitational collapse the shape 
of the turbulence induced lognormal PDF begins to be¬ 
come skewed toward the hi gh density material which 
mani fest as a power-law tail (lKlessenll2000l : ICollins et al.1 
I20I2I1 . Furthermore, the PDF was shown to be im¬ 
portant for analytic model s of star formation rates 
and initial mass functions (iPad oan &: Nordlundj 


Krumholz&^cKee 120051: iHennebelle fc Chabried 


Padqan fc Nordlund I20III : iFederrath fc KlessenI l2012i 


2002 


2008 


What can be learned from the lognormal (i.e. turbu¬ 
lence dominated) portion of the PDF? Several authors 
have suggested the turbulent sonic Mach number can 
be estimated fro m the calculation of the density/column 
density variance ([Padoan et al.1 119971 : 1 Price et all 1201111 
and the density/column densi ty skewness and kurto- 
sis, i.e. higher order moments (|Kowal fc Lazari^ 120071 : 
IBurkhart et al.ll2009[ l. 

In particular, the relationship between Ais and the 
variance of the logarithm of the 3D density distribu¬ 
tion as seen in numerical models ([Padoan et al.1 119971 : 
IPassot fc Vazauez-Semadenil 1199811 generally takes the 
form: 

( 10 ) 


'-P/Po = 


where po is the mean value of the 3D density field, b is 
a constant which depends on the driving of the turbu¬ 


lence in question with 6=1/3 for solenoidal forcing and 
6=1 for compressive driving, and cr is the standard de- 
viation of the density field normalized by its mean value 
(iNordlund fc Padoanlll99ll IFederrath et al.l l2008l 120101 : 
IMolina et al.l 1201^ 

When taking the logarithm of the normalized density 
field this relationship becomes: 

a!=Hl + b^Ml) ( 11 ) 

where s = ln(p/po) and cr^ is the standard deviation of 
the logarithm of density (not to be confused with Up/po). 

The relationship between the observable column 
density standard deviation and sonic Mach number 
([Burkhart &: LazarianI |2012[) retains the same form as 
that of the 3D density field but with a scaling constant 
A: 


= A\n{l + b^Ml) (12) 

The corresponding relation for the linear variance 
based on Equation [TS] is: 

4/Eo = + 1)^-1 (13) 

where C, = ln(I]/I]o) and E is the column density dis¬ 
tribution with Eq denoting the mean value of the column 
density distribution. 

In addition, the PDF of incremental fluctuations has 
been shown to be useful for studies of turbulence in 
the density regimes where gravity is not dom inant 
([Esquivel fc Lazarianll20lTil : iTofflemire et aLll201Ifl . The 
Tsallis statistic provides a fit for incremental PDEs and 
the fit parameters describing the width and amplitude 
are related to the physics of the gas such as the sonic 
and Alfvenic Mach numbers. 

In the high density regime (i.e., A y > 2) regime. 
where gravity dominates the PDF shape ([Schneider et al.1 
I2m4all . power law tails for m with exponent related t o 
the star formation efficiency ([Federrath fc Klessenll20I^ . 
These power law tails have been observed in n umerical 
simulations of both density (IC ollins et al.l l20I^ and col¬ 
umn density (IFederrath &: Klessenll20I3ll as well as ob- 
servatio ns of dust extinction in numerous star formin g 
regions ([Kainulainen et aLll2nill : ISchneider et al.ll2ni4all . 

In this section we focus on the column density PDFs, 
Tsallis function and PDF moments in MHD turbulence 
simulations with and without gravity for a range of mag¬ 
netic field strengths. In the case of the Enzo simulations, 
we pay particular attention to the time evolution in order 
to see how the column density PDE evolves as the cloud 
collapses. 


3.1. Column Density PDF 

Eigure|2]shows the column density PDEs of the Enzo 
simulations for different line-of-sight (LOS) orientations 
(rows across) and different time steps (in units of t//) 
indicated in different colors/linestyles, with the t = 0 
case indicating pure supersonic MHD turbulence plotted 
as a solid black line. The top row displays LOS in the x- 
direction (i.e. the direction of the mean magnetic field), 
the middle row displays the LOS in the y-direction, while 
the bottom row displays the LOS in the z-direction. 

We fit lognormals to the peaks of the distribution and 
separately fit a power law for t > 0.35t// in the column 
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-1 0 1 234 -1 0 1 234 -101234 


Figure 2. PDFs of InS/Eo (i.e. () at 512® resolution. The top row displays line of sight (LOS) in the x-direction (i.e. the direction of 
the mean magnetic field), the middle row displays the LOS in the y-direction, while the bottom row displays the LOS in the ^-direction. 
Columns show the high, mid and low Bg, respectively. Different time steps (in units of tg) are shown with the color and linestyle given in 
the legend. We overplot the density regime that we fit the power law tails as thick solid lines for t > 0.35tyy. 


density range of C = 1-3 — 1.9. This fit range avoids the 
very high density portion of the PDF, which is plagued 
by low number statistics, as well as the peak of the dis¬ 
tribution. This range is denoted by the thicker red line 
in plot. We divide the PDFs into three different time 
regimes and point out the following visual features: 

1. At t < 0.15t//, PDFs (black solid and yellow 
dashed lines) display general lognormal behavior. 
The simulation at t = 0.15t// shows a wider log 
normal than the t=0 simulation for roughly the 
same sonic Mach number. For t < O.lbtff in the 
high density range we do not observe a clear power 
law tail or the tail is very steep such that it is vi¬ 
sually indistinguishable from the lognormal 

2. At 0.35t// <t< 0.45t//, the PDFs retain the gen¬ 
eral log normal shape for log column density values 
C < 1 regardless of LOS orientation or magnetic 


field strength, however for column densities C > 1 
power law tails are clearly seen for time steps at 
t > 0.35t// 

3. At t > 0.45t//, the PDFs (purple dotted lines) 
show clear power law tails at the high column den¬ 
sity end while the PDF distribution of low column 
density material has become broadened as com¬ 
pared with earlier times. This broadening of the 
low column de nsity material has als o been observed 
in density Isee lCollins et al.l (j2012t) ). 

Examining the PDFs of Figure [2] also shows that there 
are differences in the behavior for different LOS orienta¬ 
tions and that the magnetic field strength plays a role in 
the PDF shape and onset of the power law tail. In par¬ 
ticular, turbulence with a lower magnetic field (far right 
column denoted with low in the title of the plot) ex¬ 
hibits wider PDFs at low densities and more pronounced 
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and shallower power law tails at the higher densities. 
Both the wider PDFs and shallower power law tails are 
more exaggerated when the LOS is parallel (top row) to 
the magnetic held rather then perpendicular (middle and 
bottom rows). We investigate these effect quantitatively 
in the following subsections by exploring the PDF mo¬ 
ments and Tsallis statistics as well as fitting the power 
law tail exponential, a. 


3.2. Power Law Tail 

The PDF of the collapsing matter forms a power 
law, i.e. y{p) ~ p “ , for densities above 
a critical density (iKlesse 3 tlOOl ISlvz et al.l [20?)^ 
IVazauez-Semadeni et al.ll200^. H igh resolution AMR 
simulations bv iKritsuk et al.l (|201lll measured a range of 
slopes: -1.67 at intermediate densities and -1.5 at high 
densities. They posit that the flattening of the slope is 
due to the onset of rotational support, which is backed 
up by the analysis of th e support function of that data 
bviSchmidt et al.l (I2013I1 . Recently, iFederrath fc KlessenI 
(I2oT1 ~ investigated the power law tail indices of self- 
gravitating MHD simulations and found that the high 
density tails are consistent with equivalent radial den¬ 
sity profiles, p oc r~'^ with k = 1.5 — 2.5. Observational 
constraints of the power law tail of the column density 
PDF are comparable with ranges reported by simu lations 
(|Arzoumaniaii et ^120111 : iSchneider et al.ll2014all . 

We plot the power law tail exponent a versus time 
evolution for t > 0.2t/y averaged over the three cardinal 
LOS directions with error bars indicating the standard 
deviation between different LOS in Figure O The LOS 
relative to the mean magnetic field is generally an un¬ 
known parameter for observational studies. We fit the 
slope to the range of column densities show in Figure [2] 
(i.e. InE/So = 1.3 — 1.9). These values were chosen 
to avoid both the low number statistics of the high den¬ 
sity portion of the PDF and the peak of the PDF. We 
will address a robust fit for the power law tails based on 
maximum-likelihood htting in a future work. 

We include the published values for the power law tail 
slopes and approximate number of young stellar objects 
(YSOs) for five molecu l ar cloud s from Herschel data from 
ISchneider et al.l (|2013l . I2014b[l : NGC3603, which is an 
active low-mass star forming region the Auriga cloud, as 
well as the Orion B, Maddalena, and Aquila clouds. We 
discuss these clouds and their relation to the simulations 
in the discussion section. We omit the PDF estimates 
from the Carina cloud as the dynamics of Carina may 
include significant massive stellar feedback which is not 
treated in our simulatio ns. This comparison is comple - 
mentary to a study by iKainulainen fc Hennind (|2014r) 
who showed that a proxy for the 3D density PDF could 
be related to the number of YSOs in several different 
clouds. 

The power law tail index is steepest for earlier time 
steps and generally becomes increasingly shallow as the 
cloud proceeds to collapse. This is expected from past 
studies of the column density powe r law tail index and 
is also observed in the density PDF (iKritsu k et al.|| 2011l : 
iCollins et al.ll20l3 : IFederrath fc Klessenll2012i 120131 1. fir 
addition to the shallowing of the power law index with 
time. Figure |3] shows that there is a strong dependency 
on the magnetic field. Simulations with higher field 
strengths show steeper slopes for all times steps due to 


the suppression of density enhancements by the magnetic 
field. The low magnetic field run shows shallower values 
of a across the time evolution parameter space. 

The low magnetization simulation has values of a rang¬ 
ing from -2.5 to -0.7, which falls out of the bounds set by 
the observations at later times. This simulation has an 
Alfvenic Mach number k, 40, which may be two high for 
realistic clouds. This may suggest super-Alfvenic clouds 
have Alfvenic Mach numbers in the range of 6-12, which 
is representative of our high and mid magnetic field sim¬ 
ulations a 2 ff_Jias_been_aj£ued_for in a number of past 
works dPadoan fc Nordhmdl 119991: iLunttila et al.l 120081: 

Burkhart et al.l 120091 : iCrutcher et al.ll2009HCollins et al.l 

2nilL 


3.3. Column Density PDF Variance 

We investigate the PDF variance as a function of time 
evolution in Figure 0] for the Enzo and Codunov simu¬ 
lations. We investigate several different methods of cal¬ 
culation of the variance in Figure |4j including directly 
calculating the variance of the column density distribu¬ 
tion (cl/^odir’ directly calculating the 

variance of the natural logarithm of the column density 
distribution center panel). 

For a distribution E = = l...A^ the mean value 

is defined as: Eq = (^*)- We define the direct 

calculation of the variance as: 


O’s/Sodir — 




(14) 


The calculation of the variance directly from the data 
does not assume that the PDF follows a particular model 
(i.e. that it is a lognormal). We also calculate directly 
the variance of natural logarithm of the column density 
distribution as: 


1 

E (C* - Co)' (15) 

i=l 

Figured! bottom panel, shows the variance calculated 
by fitting a Caussian to the peak of the distribution and 
measuring the variance from the fit (ct^^^j). 

The column density variance of the MHD simulations 
at t = 0 depends primarily on the sonic Mach num¬ 
ber (iBurkhart fc Lazarianl[2M^ regardless of the calcu¬ 
lation method. The distribution of gas with larger sonic 
Mach number shows increasing variance for both direct 
calculation of the variance and the variance calculated 
from a Caussian fit. When gravity begins to alter the 
PDF from lognormal, i.e. at t > 0, the variance 

and dir) 1^® Enzo simulations (which have sonic 
Mach numbers of ~ 10) dramatically increases past the 
expectations of MHD turbulence even with Mach num¬ 
bers as high as 20. This effect manifests as several orders 
of magnitude larger difference in the directly calculated 
variance (top panel) of the column density distribution. 
The natural logarithm of the column density distribution 
(middle panel) also shows a linear increase in variance 
as the cloud evolves with gravity. The variance as ob¬ 
tained from a Caussian fit remains generally flat within 
the error bars across the time evolution parameter space 
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Figure 3. The power law tail index (o) of the 512® resolution Enzo simulations, averaged ove r all LOS, versus time evolu tion. Error bars 
denote the standard deviation between different LOS. Lines indicate observed powerlaws from ISchneider et al.l 112013 1 and lSchneider et al.l 
Il2014bll 


up until about t = QAIq for the low and mid magnetic 
field runs and for the high field run the variance does 
not significantly change in the time parameter space in¬ 
vestigated. This is because the Gaussian fit is applied 
to the peak of the distribution and is sensitive primarily 
to the lower density material dominated by turbulence 
and not the power law tail high density portion of the 
PDF. However, there is a slight trend toward a wider 
lognormal at the low density end of the PDF (which can 
be visually seen in Figure [5]). This suggests that fitting 
a Gaussian to the peak of the distribution can be used 
to dissect the turbulence dominated portions of the gas 
while the fraction of gas collapsing to form proto-stars 
may be probed with the formation of the power law tail 
index (see Figure O. 

Additionally the variance of collapsing column density 
material shows a dependency on the magnetic field. For 
all three methods of variance calculation, the simulation 
with lower field strength shows higher values of variance 
as compared to the high/mid cases which are generally 
not distinguishable within the error bars. This effect 
is more pronounced and increasingly significant across 
different sight lines as the time evolution increases. As 
the magnetic field increases so does the suppression of 
high density clumps, which causes the column density 
variance to be smaller compared with simulations with 
lower field strengths. The dependency of the column 
density variance on cloud magnetization in gravoturbu- 
lence simulations is somewhat surprising as the variance 
has been shown to be weakly dependent magnetic field 
strength in other works that focus ed only on supersonic 
MHD simulations without gravity (jBurkhart et al.ir201(1l : 


iMolina et al.ll201^ . These results suggest the magnetic 
field plays a role in the global support of the cloud against 
gravity. 

Gomparison wit h published dust ex t inctio n column 
density maps from iKainulainen &: T^ (I2013D (Table 1 
of that work) shows that clouds have values of directly 
calculated variance ranging from ~ 0.64. 

Taking the logarithm of these numbers for ease of com¬ 
parison with the top panel of Figure 21 these values are 
loga^yj,^ = —0.6 to —0.19, which suggests the cloud vari¬ 
ance can not be attributed to MHD turbulence alone, and 
that gravity must be active up to t = 0.5t// as compared 
with ou r high or mid magnetic field simulations. We over 
plot the IKainulainen fc Tad (|2013ll cloud variance ranges 
with straight dashed lines in the top panel of Figure 2) 
We compare the self-gravitating snapshots from the 
Enzo simulations to the column density variance - sonic 
Mach number relation given in Equation (1131) in Eigure[5j 
Our general conclusions indicate that two regimes exist 
which are defined by the importance of self-gravity: 

1. The column density variance - sonic Mach num¬ 
ber relation without gravity: the Godunov (black 
squares) and t = 0 Enzo simulations (black symbols 
denoted in the legend) follo w closely the predic¬ 
tion o f equation 231 given by iBurkhart &: LazarianI 

(EoHI). 

2. The column density variance - sonic Mach number 
relation with gravity: once gravity becomes impor¬ 
tant the PDF variance no longer tracks the behav¬ 
ior of the sonic Mach number. The more evolved 
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Figure 4. Different methods of calculation of the PDF vari¬ 
ance versus time. The Godunov MHD turbulence simulations 
are plotted with green symbols at t=0. Top row: directly calcu¬ 
lated variance of the linear column density distribution. Straight 
dashed lines acro s s rep resent the range of values found in the 
IKainulainen et ^ 1I2009I~) survey of IRDC PDFs. Middle row: di¬ 
rectly calculated variance of the natural logarithm of the column 
density distribution. Bottom row: Gaussian fitted variance of the 
column density distribution. Error bars are calculated as the stan¬ 
dard deviation between three different sight lines. 
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the cloud’s collapse, the higher the variance even 
for the same sonic Mach number. 

Figure [5] highlights the importance of the magnetic 
state of the gas in the evolution of the PDF of GMCs. 
The variance of the Enzo simulations depends strongly 
on the global strength of the magnetic field. The lower 
the value of the magnetic field, the higher the measured 
variance. This is because simulations with stronger field 
suppress the formation of dense cores as a function of 
time. 

Figure [S] also compares the Enzo and Godunov col¬ 
umn density variance - sonic Mach number rela tion with 
that of the GMG s fromiKain ulainen et al.l (1200911 as taken 
from Table 3 of IKainulainen fc TanI (j2013fl . The GMC 
variance values are all larger than the Godunov simula¬ 
tions and follow the upward evolution of the Enzo simula¬ 
tions along the variance axis. A direct comparison of the 
range of variance values with clouds that have Mg ~ 10 
suggests evolution time scales of f = 0.25 — 0.45f// over 
a range of magnetic field values. 

3.4. Column Density PDF Higher Order Moments 

In addition to the variance, the higher order moments 
of the PDF have also been shown in a number of works 
to be sensitive diagnostics of the sonic Mach number 
of a turbulent density fie ld (|Kowal fc Lazari^ 120071 : 
iBurkhart et al.l[2009L 1201011 . However these earlier works 
investigated the higher order moments of the linear col¬ 
umn density distribution (i.e. E/Sq). As shown in the 
previous subsections, the power law tail of the column 
density distribution is a sensitive diagnostic of the evolu¬ 
tionary state of a collapsing cloud. A additional sensitive 
diagnostic that could be complimentary to the power law 
tails seen in the natural logarithm of the column density 
distribution ((^) are the higher order moments precisely 
because they characterize deviations from Gaussianity. 

Skewness and kurtosis are defined by the third and 
fourth-order statistical moment. Skewness is defined as: 



If a distribution is Gaussian, the skewness is zero. Kurto¬ 
sis is a measure of whether a quantity has a distribution 
that is peaked or flattened compared to a normal Gaus¬ 
sian distribution. If a data set has positive kurtosis then 
it will have a distinct sharp peak near the mean and have 
elongated tails. If a data set has negative kurtosis then 
it will be flat at the mean. Kurtosis is defined as: 



We plot the higher order moments of the natural log¬ 
arithm of the column density distribution in Figure |6l 
The expectations for MHD turbulence with no gravity 
(i.e. our t = 0 cases) are values of skewness and kur¬ 
tosis bounded around zero as the lognormal distribution 
should have relatively small skewness and kurtosis for 
given sonic Mach numbei0. An increase in sonic Mach 
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This is not the case for the linear distribution of density. For 
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Figure 5. Godunov simulations and Enzo simulations. Pink, green, orange and red colors represent increaseing 
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Figure 6. Higher order moments vs. time evolution of the natural 
logarithm of the column density distribution (i^). Skewness (5) is 
shown in the top panel and kurtosis (K) is shown in the bottom 
panel. We use the same color scheme as Figure [4| Later times 
show higher skewness and kurtosis but unlike the variance, little 
dependency magnetic field variation is found. 


example. IBurkhart et aP (120091) showed that the skewness and kur¬ 
tosis are sensitive to sonic Mach number for linear density and 
column density as shocks increase tails toward the high density 
portion of the PDF 


number affects the width of the lognormal PDF and gen¬ 
erally not the tails or peaks. Indeed, at t < 0.2tff the 
values of skewness are bounded between 0.4 and -0.4 and 
the values for kurtosis are bounded between -0.6 and 0.8. 
For t > 0.2tff the values of skewness and kurtosis are 
generally positive and increase with time evolution due 
to the formation of the power law tail that creates more 
peaked distributions and skews the PDF tail toward the 
high density end. In general, we do not see a systematic 
sensitivity to the magnetization of the simulations, how¬ 
ever the low and mid simulations generally show higher 
values of skewness and kurtosis than the simulation with 
high value of magnetic field. 

Figured] suggests that calculation of the skewness and 
kurtosis of the natural logarithm of the column density 
distribution (i.e. C) can be complimentary to fitting the 
power law tails to the high density material. We test this 
idea by plotting the power law tail index a vs. the skew¬ 
ness and kurtosis of C in Figure |7l Both the higher order 
moment and the power law tail index increase mono- 
tonically with time evolution and hence can be seen to 
increase together. Although the power law tails have de¬ 
pendency on the magnetization of the gas (see Figure [3]), 
the higher order moments of C are somewhat insensitive 
to field strength. This is clear in that the lowest magnetic 
field run (blue data) proceeds out to much higher values 
of a for the same spread in values of skewness/kurtosis. 
This implies that the degeneracy can be broken and that 
observers who measure higher order moments and power 
law tails out to larger values might constrain both the 
energetic importance of the magnetic field as well as the 
time evolution of clouds. 

The many avenues of exploring PDFs, including fitting 
a Gaussian to the peak and a power law to the high den¬ 
sity tail region, as well as direct calculation of the PDF 
moments of the entire distribution provide researchers 
with a tool kit for separating out the contributions of 
MHD turbulence and gravity in the structure formation 
and evolution of clouds. In the next sub-section we also 
suggest spatial incremental PDFs to be of use toward 
this purpose. 
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3.5. Column Density Incremental PDFs: the Tsallis 
Statistic 

In studying PDFs we study the column density field 
as it is given from observations. The density may be 
influenced by multiple processes that act differently on 
different scales. Thus it is advantageous to use a measure 
which can differentiate the properties of turbulence at 
different scales, e.g., at the scales of energy injection, 
inertial range interval, energy dissipation scale and scales 
where gravitational interactions become important. To 
do this, it is advantageous to use incremental measures 
that are describe the increments of densities over the 
scale r, namely, pix + r) — p{x). Such incremental PDFs 


may be studied by different means and fitting those to 
the Tsallis distribution is a way to do the job. 

The Tsallis distribution is often used in the context 
of non-exte nsive statistica l dynamics. It was originally 
derived (see iTsaili^ (jI988[) 1 from an entropy generaliza¬ 
tion to extend the traditional Boltzmann-Gibbs statistics 
to multi-fractal systems (such as the ISM). The Tsal¬ 
lis distribution has since been applied to problems over 
a range of applications as the distribution can describe 
PDFs with tails that are not exponentially bounded. 

The Tsallis function of an arbitrary incremental PDF 
(A/(r)) has the form: 
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Rq = a 


l + (<7-l) 


A/(r)2 






(18) 


where A f{r) denotes the normalized PDF of incremen¬ 
tal differences (with spatial separation r), i.e. 


A/(r) = (/(x,r) - (/(x,r))^)/cr/, 


(19) 


where cr/ the dispersion of the increments and (...)x 
denotes spatial averaging over a shell of size r = |r|. 
The normalization used is such that the PDF has mean 
value at zero and a standard deviation of unity. We 
denote the column density increments as: /(x, r) = 
S(x r) — S(x), A normalized histogram of the incre¬ 
mental maps for a given lag results in our incremental 
PDF which is then fit with the Tsallis function. 

The other parameters in Eq. (HI are as follows: q 
is the so called “entropic index” or “non-extensivity pa¬ 
rameter” which is related to the size of the tail of the dis¬ 
tribution; w is a measure of the with of the PDF related 
to the width of the distribution; and a, the amplitude. 
By varying the parameter q in the Tsallis distribution it 
is possible to obtain distributions that range from Gaus¬ 
sian at g —>■ 1 to “peaky” distributions with large tails. 
The parameter q is closely related to the kurtosis (fourth 
order one-point moment) of the PDF, and similarly the 
parameter w is related to the variance of the PDF. 

The Tsallis distribution reduces to the classical Boltz¬ 
mann - Gibbs (Gaussian) distribution in the limit of 
(7 —>■ 1. However, for the present purpose we use it 
as an empiric function that fits well the properties of 
MHD turbulenc e as was shown i n lEsquivel &: LazarianI 
(see also iTofflemire et "al] (|201ll) '). We also note 
that the Tsallis function was successfully used to char¬ 
acterize the magnetic field of the sol a r wind in a series 
of pap e rs bv iBurlaga fc Viiiasl (I2004D : iBurlaga fc Vinai 


(2004); 

(2005) 


Burlaga fc -VifiasI (120051): Burlaga fc F.-VinAs 

Burlaga et ahl (I200f)l 1200711: iBurlagal (1200.111 


lEsauivel fc LazarianI (I2010D and fTofflemire et akl ( 201 ll) 
used the Tsallis statistics to describe the spatial variation 
in PDFs of turbulent density, column density, velocity, 
a nd magnetization of MHD sim ulations without gravity. 

lEsq^uivel fc LazarianI (|2010D and ITofflemire et al.1 
(|2011l) used the Tsallis statistics to describe the spatial 
variation in PDFs of turbulent density, column density, 
velocity, and magnetization of MHD simulations with¬ 
out gravity. Both efforts found that Tsallis provided 
adequate fits to their incremental PDFs and gave insight 
into statistics of MHD turbulence. Our present study 
is the first attempt to apply the Tsallis statistics to the 
density field obtained with the simulations that include 
self-gravity. 

F or non self-gravitating MHD simulati ons it was shown 
in lEsauivel fc LazarianI (j20ir)H and ITofflemire et al.1 
(|201ll) that the Tsallis fit parameters a, q and w of the 
column density distribution were sensitive to both the 
sonic and Alfven Mach numbers. Higher sonic Mach 
numbers and higher magnetic field strengths produced 
incremental PDFs with higher values of width (re), am¬ 
plitude (a) and kurtotic index (q). In this subsection we 
investigate the first use of the Tsallis function to simula¬ 
tions of self-gravitating MHD turbulence. 

Figure H] shows the three Tsallis parameters, a, w and 
q vs. spatial lag (in pixels) for the Enzo simulations. 


Error bars are plotted by taking the standard deviation 
of values of the fit parameters with different LOS rel¬ 
ative to the mean magnetic field. High, mid, and low 
values of magnetic field are presented for comparison in 
the columns going from left to right, respectively. In 
simulations with and without self-gravity, signs of the 
dissipation scales can be seen at small lag increments. 

The solid black line represents the Enzo simulation 
with pure supersonic MHD tur bulence. Gompar i son o f 
these values with Figure 8 from ITofflemire et al.1 (j2011f ) 
show very good agreement with other supersonic MHD 
simulations. Once self-gravity begins to create regions 
of over-density, all three Tsallis parameters dramatically 
increase well past the expectations for supersonic turbu¬ 
lence. 

The values of the amplitude and width of the Tsallis 
column density PDFs show dependencies on the mag¬ 
netization of the gas. For the low and mid magnetic 
field simulations, the incremental PDFs are wider and 
have higher amplitudes compared with the high mag- 
ne tization case. This is in contrast to the findings 
of ITofflemire et al.1 (j201ir i , which found that a higher 
field increases the Tsallis parameters in MHD turbulence 
without gravity. However, when the material collapses 
the magnetic field acts to suppress overdense regions 
from forming which in turn creates incremental PDFs 
with lower values of width and amplitude. Within the 
error bars given by the LOS orientation the effect of the 
magnetic field is not distinguishable between the mid and 
low magnetic field simulations but is nearly an order of 
magnitude different comparing these two regimes with 
the high magnetic field simulation at later times (i.e. 

The kurtotic index q shows little clear variation with 
changing magnetic field but does show a dependency on 
the time evolution of the column density. As gravity 
acts to create the dense clumps seen in Figure 1, the 
kurtotic index q is seen to increase past the turbulence 
only snapshot (black line). The increase in q is mono¬ 
tonic with increasing time step. The fact that q does 
not depend strongly on magnetic field (or sonic Mach 
number, see (ITofflemire et al.ll2011h l but rather only on 
the collapse evolution suggests that this parameter might 
be extremely useful in breaking the degeneracy between 
gravity, compressibility and magnetization in the PDF 
statistics of column density. We discuss this further in 
section O 

4. COLUMN DENSITY POWER SPECTRA 

Complementary to PDFs, an essential tool for tur¬ 
bulence studies is the spatial density and velocity 
power spectrum. The turbulence energy transfer pro¬ 
cess can be studied by examining the Fourier power 
spectrum, and the sources and sinks of energy, in¬ 
cluding the injection and dissipation scales, can be 
identified. The power spectra of density and veloc¬ 
ity (and their variants such as the structure func¬ 
tion and delta variance) have been suggested by sev¬ 
eral authors to provide information on the spatial 
and kinematic scaling of turbulence, sonic Mach num¬ 
ber a n d injection/dissi p ation sc ales (iK owal fc LazarianI 
20071: iBurkhart et al.l I^OlOl: lOssenkoof et al.l 120011 : 
Collins et al.ll2012l : iFederrath fc KlessenI I2013D . In this 
section we explore the use of the spatial power spectrum 
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Figure 8. Tsallis fit parameters vs. incremental lag. The top row shows the logarithm of the amplitude parameter (Log a), the middle 
row shows the logarithm of the width parameter (Log w) and the bottom row shows the kurtotic parameter (g). Columns show high, mid 
and low magnetic field strength simulations, respectively. The color scheme follows that of Figure [2 except that we include two additional 
models: purple solid line with t = 0.55tyy and yellow dashed and dotted line with t = 0.25£//. 
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of column density maps of self-gravitating MHD turbu¬ 
lence for determining the dynamical evolution of clouds 
undergoing collapse and investigate if the power spec¬ 
trum might be a complimentary tool to the PDF. We 
also investigate the origin of the changes in the column 
density power spectrum as a function of time. 

The Fourier transform of the two point autocorrelation 
function (i.e. the spatial power spectrum) provides in¬ 
formation on the properties of the turbulence cascade, 
including the injection scales and dissipation scales. The 
power spectrum is defined as: 

P(k)= ^ F(k)-F*ik) (20) 

]<.—const. 

where k is the wavenumber and F(k) is the Fourier trans¬ 
form of the field under study, which for our purposes is 
the synthetic column density maps. 

The one-dimensional energy spectrum E{k) is re¬ 
lated to the measured power spectrum by E{k)dk oc 
P{k)dk^, where D is the dimensionality. For incom¬ 
pressible turbulence, the Kolmogorov power spectrum 
scaling (jKolmogorovIllO^ in three dimensions (3D) is 
P{k) 3 D oc and the energy spectrum scales as 

E{k) oc k~^/^. For the same E{k), P{k) 2 D oc 
and in ID P{k)iD oc k~^/^. 

Although the Kolmogorov slope is for incompressible 
unmag netized fluids, the analysis of iGoldreich fc Sridharl 
(11995(1 showed that the energy spectrum scaling of 
the incompressible cascade perpendicular to the mean 
magnetic field also retain s the —5/3 slope. The 
Goldreich &: Sridliail (1199^ an alysis was exteiided in 
Lazarian fc Vishnia (j 1999( 1 and iGho fc Vishnia (|200flf l 
to include the concept of the cascade relative to the 
local magnetic field to obtain the correct scaling re¬ 
lations. Later studies confirme d the —5/3 slope with 
highe r resolution simu lations (|Beresnvak fc LazarianI 
(20091 : (Beresnv^ (2012(1 . The actual spectrum of 
MHD turbulence is anisotropic, and scale-dependent 
anisotropy in the system of reference is connected with 
the local magnetic field (see the discu ssion of the con¬ 
cept o f the local magnetic field, e.g. see [Cho &: Lazari^ 
((2003(1 1. The statistics of fluctuations given by the spec¬ 
tral slope is different when measurement are made par¬ 
allel and perpendicular to the local magnetic field. How¬ 
ever, in this paper we deal with LOS observations and in 
this case it is di fficult to measure such anisotropy in in¬ 
tensity maps (see [Burkhart et al.[ (|2014[1 for a discussion 
of how this can be done in velocity centroid maps). For 
our purposes in this work we can use the —5/3 reference 
slope as we deal with LOS observations of super-Alfvenic 
turbulence and not the spectrum as measured in the local 
frame of reference to the magnetic field. 

In the current paper we are interested in the behavior 
of the density/column density spectral slope in compress¬ 
ible self-gravitating turbulence. In the presence of su¬ 
personic turbulence, such as exists in GMCs, the density 
spectral slope is shallower than the relations discussed 
above d ue to shocks creating small scale enhancements o f 
density ((Beresnvak et al.[r2005KKowal fc LazariaD([200^ . 

[Burkhart et al.[ (|2010(l plotted the power spectral slope 
of column density maps versus sonic Mach number and 
found that the slope of the power spectrum of ideal MHD 


turbulence is increasingly shallow as the Mach number 
increase. However they found that the power spectral 
slope begins to saturate toward -2 (as compared with 
the 2D slope of -8/3) for very high sonic Mach number, 
regardless of Alfvenic Mach number. 

Furthermore, the inclusion of gravity in a magne¬ 
tized turbulent media changes the behavior of the 
density spectral slope dramatically. As gravity fur¬ 
ther enhances overdense regions, the spectrum be¬ 
comes increasingl y shallow as more ma t erial collects 
on small scales ([Ossenk opf et al.l [2001[ : [Collins et al.[ 
[2012[: (Federrath fc KlessCT (201^ . In some cases, 
self-gravitating supersonic turbulence produces density 
structure that drive the spectral slope toward positive 
values. This is in contrast to non-self-gravitating turbu¬ 
lence where the power is dominated by large scale struc¬ 
tures and the power on the smaller scales is decreasing. 

In this section, we investigate the evolution of the ID 
column density power spectrum (which we denote as 
P{k)) as one transitions from gas dynamics dominated 
by supersonic MHD turbulence to self-gravitating. We 
also propose a model for the power spectrum of a self- 
gravitating fluid based on our findings. For the AMR 
data, the analysis was done on a coarse-grained model, 
where the refinement was restricted via volume-weighted 
average to the coarsest level. This is due to the sparse 
sampling of data at these higher wave numbers, which 
would lead to unphysical suppression of power at these 
wave numbers. 

We plot the ID power spectrum averaged over three 
lines of sight as a function of wavenumber in Figure IHl for 
the Enzo simulations. We overplot a solid black line in 
our inertial range with slope of -5/3 from log fc=l.l to 
log fc=1.4 for reference. Solid colored lines are overplot¬ 
ted curves of Equation (next subsection) with least 
square fit parameters given in Table [2] for each model. 
We find that there is no significant variation of the col¬ 
umn density spectra with projection direc t ion, s imilar to 
past studies such as [Federrath &: Klessed (|2013() . 

Similar to the PDFs, we find that there exists three 
distinct stages of evolution in the power spectrum of col¬ 
lapsing column density images: 

1. At t < 0.15tff (henceforth termed “early”), the 
cloud is in a purely turbulent regime, and hence 
the power spectrum exhibits behavior of supersonic 
turbulence, i.e. negative valued slopes which are 
shallower than the -5/3 slope. 

2. At 0.15t// < t < 0.35t// (henceforth termed “in¬ 
termediate”): As the timestep increase and grav¬ 
ity begins to dominate the small scales (large fc), 
the slope becomes increasingly shallow and at the 
largest timesteps shown, is positive in the inertial 
range. This is in agreement with past studies of 
the powe r spectrum of collapsing supersonic tur¬ 
bulen ce (jGollins et al.[ \201‘A [Federrath fc KlesserJ 
[2013( 1. The turn over time step from negative to 
positive occurs around t ~ 0.25t// but is also de¬ 
pendent on the strength of the magnetic field. 

3. At t > 0.35tff (henceforth termed “advanced”), 
gravity dominates the power spectrum at large k 
values and the slopes in both the inertial rage and 
at larger k are positive valued. 
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Figure 9. ID Power spectra as a function of wavenumber, av¬ 
eraged over three lines of sight, for the Enzo simulations at reso¬ 
lution 512^. Top to bottom plots represent high to low magnetic 
field runs, r espe ctively. Solid colored lines are overplotted curves 
of Equation 1211 with least square fit parameters given in Table [2] 
for each model. The color scheme is identical to Figure 7. 


It is interesting that the column density power spec¬ 
trum (which is an observable quantity) exhibits different 
behavior based on time evolution. This implies that this 
tool maybe used on observations to not only determine 
the properties of turbulence in clouds but also the grav¬ 
itational state of cores. We develop a functional fit to 
determine the turnover scale at intermediate times tran¬ 
sitioning from turbulence dominated to gravitationally 
dominated in the next subsection. 

We plot the power law slopes (denoted as /3i) of the 
spectrum (as shown in Figure[9]) versus time for the Enzo 
and Godunov (green symbols) simulations in Figure fTOl 
Error bars are calculated by taking parallel and perpen¬ 
dicular sight lines relative to the mean magnetic field. 
For each sonic Mach number of the Godunov simulations, 
two separate Alfvenic Mach numbers are plotted. 

The Enzo simulations at t = 0 show slopes that are 
fully consistent with the Godunov MHD simulations for 
the same sonic Mach number. For t > 0 the slopes of 
the Enzo simulations begin to become increasingly shal¬ 
low as compared with the purely turbulent scenarios, as 
is evident in Figure [9l However this is not significant 
within the error bars until t > 0.2tff. The most sig¬ 
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Figure 10. Column density power spectral slopes as a function of 
time for the Godunov and Enzo simulations with varying values of 
magnetic field. The color-symbols used is similar to Figure 4. Er¬ 
ror bars are created from the standard deviation of different sight 
lines. The MHD simulations show the expected slopes for subsonic 
values around —1.6 and becoming increasingly shallow for su¬ 
personic turbulence. The Enzo simulation at t=0 are within the 
expected slope range of the Godunov/ENO simulations however, 
once gravity is turned on the values of the slope increase past the 
purely supersonic turbulence cases and eventually become positive. 


nificant gains in the value of the slope are made in the 
case of the low magnetic field. At t > 0.25t// this slope 
becomes positive and remains so as the cloud continues 
to evolve. For the high and mid magnetic field cases the 
column density power spectral slope becomes positive at 
fa t > 0.45Gf. In both this work and previous works 
such as iFederrath fc KlessenI (|2013ll . positive values for 
the column density spectral slope are observed at evolved 
collapse stages. 

4.1. Modeling the Power Spectrum of Self-gravitating 
Turbulence 

Figure |9] shows that the power spectrum of the gas 
column density map reveals different behavior for super¬ 
sonic MHD turbulence and supersonic MHD turbulence 
undergoing collapse. Namely that supersonic turbulence 
displays a negative power spectral slope while in the col¬ 
lapsing state the power spectrum becomes increasingly 
shallow past the expectations of supersonic turbulence 
and even becomes positive valued. Thi s has also been 
confirmed in other studies. iCollins et al.l (e.g. I2012h and 
IFederrath fc KlessenI (1201311 . 

In this subsection we address the turnover scale that 
occurs as the spectrum transitions from turbulence domi¬ 
nated to gravitational dominated in order to understand 
the nature of this transition and to provide observers 
with additional methods to disentangle the dynamics of 
turbulence from cloud collapse. 

We propose a functional fit to the power spectrum in 
the form: 
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P{k) = Aik^^exp{—k/kc) (21) 

This is characterized by an amplitude Ai and power law 
behavior k^^ which dominates the small k behavior (i.e. 
turbulence dominated). The scale kc characterizes the 
turn over from turbulence dominated to self-gravity dom¬ 
inated which is observed in the power spectrum at an in¬ 
termediate time evolution (see Figure |9|). The advantage 
of the exponential form is that it can roughly describe all 
three stages of evolution discussed in the previous section 
and hence could be used by observers who do not know 
apriori the time scale. At t=0 “early” timesteps the ex¬ 
ponential can be negative to describe the dissipation of 
energy (a form used in p revious works on MHD turbu¬ 
lence, e.g. iLazariaiil (j20n4[ l ). At intermediate times it can 
describe the turnover scale seen in Figure IHl Finally at 
“advanced” times Equation [5T] reduces again to a power 
law form as kc becomes large. We fit the function from 
fc=13 to avoid the scales dominated by the turbulence 
driving. 

Figured overplots the fits of the three parameters (Ai, 
(32, and kc ) given in Equation!^ as thick solid lines. We 
present the fitted parameters in Table [5] Eor the power 
spectra with clear signatures of turbulence, i.e. power 
law behavior in the inertial range and decreasing power 
in the dissipation range (large fc), we only overplot the 
contribution of the power law in Equation 1211 We note 
this in Table [5] under the column for kc with N/A. Table 
[2] also lists the values of the slope /3i which is calculated 
from a linear fit in the inertial range from fc=13 to fc=24 
and plotted in Figure [TUI 

Figure [9] shows that Equation reasonably models 
the turn over scale at intermediate time steps and the 
gravitationally dominated high k scales (i.e. small spa¬ 
tial scales). The fits are more robust for the high and 
mid magnetic field cases, whereas the low magnetic field 
simulation (bottom panel) is not as well fit for the mid 
k ranges and hence for this simulation the values of kc 
are more erratic. In general, the 132 values from the fit of 
Equation [UT] listed in Table [U] agree within the error bar 
with the pure power law fit (i.e. /3i) for early timesteps. 

Table [U] shows that the turnover scale kc generally 
increases with increasing timestep. This is because at 
t > 0.5tff the spectrum at the large k scales has com¬ 
pletely transitioned from having a negative slope (tur¬ 
bulence dominated) to having a positive slope (gravity 
dominated). As kc increases the exponential term influ¬ 
ences Equation [UT] less until again only the power law 
term dominates again. At “advanced” time steps the 
values of kc are equivalent to infinity since they extend 
beyond the range of the power spectrum being plotted 
(i.e. past fc=256). Thus the kc scale of interest occurs at 
intermediate time steps where gravity is just starting to 
dominate the power spectrum. This is kc ~ 100—200 and 
occurs at timesteps from t = 0.25tff to t = O.dStg, with 
the difference being attributed to the influence of the 
magnetic field. Using the s c aling for these simulations 
as reported in iCollins et al.l (j2012f l. this corresponds to 
a length scale of L = 0.005 — 0.2pc. Attention to Eigure 
[U] shows that these timesteps are also where the power 
law tails begin to form in the PDEs. Visual inspection of 
the column density of these snapshots reveals that this is 
the moment when the first cores begin to form from the 
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Figure 11. Reconstructing self-gravitating spectra from turbu¬ 
lent spectra: Power spectrum as a function of wavenumber, aver¬ 
aged over three lines of sight, for the Enzo simulations at resolution 
512®. We overplot (thick solid lines) the power spectrum of the 
Enzo simulation with the to -|-max(t) column density map, i.e. the 
t=0 map that has one pixel modified such that it has amplitude 
equal to the peak value of the given snapshot at t > 0. 

turbulent cloud. Thus, the turnover scale kc can be used 
as an additional diagnostic for star formation evolution 
in a similar manner to the PDEs and the power spectral 
slope. 

4.2. Mimicking Self-gravity in the Turbulence Power 
Spectrum 

The timestep of t ~ 0.25t//, where gravity begin to 
dominate the power spectrum and PDEs of turbulent 
magnetized clouds, is the moment when very high den¬ 
sity cores begin to form in our simulations. These cores 
exist on the smallest scales (i.e. a few pixels), yet influ¬ 
ence the global turbulence statistics of the cloud. Thus 
in order to understand the affects of the small scales on 
the large scales a natural question arrises: can the global 
properties of gravity be replicated with a single (or mul¬ 
tiple) point function which replicates a core? 

We test if it is possible to replicate the power spec- 
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Table 2 

512^ Enzo power spectra fit parameters. Bext denotes the magnetic field regime (see Table 1); t denotes the time step in units of the free 
fall time; /3i is the power spectral slope fitting from a standard power law fit; Columns 4-6 (Ai, /32, ^c) show the fit parameters to Equation 
1211 Lc shows the length scale corresponding to kc assuming a 4 Mpc box Lc = 4pc * 27r/kc; (Emax) shows the ma ximu m column density, 
averaged over all three axes; PtQ+ 7 nax{t) shows the slope of the turbulence with delta function, to+n^ax(t) fSection l4.2ll : /3rec shows the fit 
to the truncated power spectrum (Section oil : Sc shows the cutoff column density in code units. We can convert the code densities to 
physical units by multiplicative scaling factor of 1.4x10^^ cm“^ assuming a cloud of mean density 1000 cm“^ and size of 4.6 pc. 


Bext 

t 


Ai 

^2 

kc 

Lc 

(Exnax) 

l^tc\-\-Tnaxit) 

/^rec 

Sc 

High 


-0.34±0.13 

0.12 

-0.5±0.03 

N/A 

N/A 

4.6 

-0.5 ±0.03 

-0.6 

4.40826 


omtff 

-0.42±0.08 

0.14 

-0.55±0.1 

N/A 

N/A 

4.8 

-0.54 ±0.0 

-0.8 

4.35521 


O.lbtff 

-0.38±0.06 

0.14 

-0.45±0.12 

N/A 

N/A 

9.9 

-0.53 ±0.01 

-0.5 

5.57656 


Q.25tff 

-0.31±0.06 

0.49 

-0.85±0.15 

120 

0.241 

65.7 

0.39 ±0.04 

-0.4 

5.41682 


O.Sbtff 

-0.36±0.01 

0.07 

-0.12±0.05 

142 

0.204 

139.7 

0.78 ±0.02 

-0.8 

5.24604 


0A5tff 

O.litO.ll 

0.01 

0.6±0.04 

700 

0.041 

176 

0.82 ±0.01 

-0.2 

5.46864 


O.bbtff 

0.31±0.07 

0.01 

0.78±0.07 

1250 

0.023 

261 

0.86 ±0.0134 

-0.1 

5.75271 


O.mtff 

0.41±0.07 

0.02 

0.76±0.06 

1452 

0.020 

290 

0.87 ±0.0 

-0.2 

5.60185 

Mid 

Otff 

-0.7±0.05 

0.12 

-0.61±0.03 

N/A 

N/A 

4.8 

-0.61 ±0.03 

-1.0 

4.58494 


O.Obtff 

-0.61±0.11 

0.09 

-0.5±0.05 

N/A 

N/A 

4.6 

-0.61 ±0.0 

-0.6 

4.10902 


0.15tff 

-0.61±0.04 

0.16 

-0.59±0.08 

N/A 

N/A 

7.4 

-0.6 ±0.0 

-0.5 

6.79376 


0.25tff 

-0.46±0.07 

0.62 

-0.95±0.03 

N/A 

N/A 

13.6 

-0.75 ±0.03 

-0.9 

10.1390 


0.35tff 

-0.27±0.0 

0.38 

-0.71±0.04 

106 

0.273 

80.3 

0.63 ±0.03 

-1.1 

5.54317 


0A5tff 

O.OlitO.16 

0.03 

0.26±0.09 

172 

0.168 

156 

0.82 ±0.01 

-0.5 

4.83240 


O.bbtff 

O.SiO.O 

0.01 

0.98±0.17 

1051 

0.028 

434 

0.88 ±0.0 

0.3 

4.60489 


Omtff 

0.68±0.09 

0.02 

0.91±0.33 

565 

0.051 

479 

0.89 ±0.0 

0.1 

4.05968 

Low 


-0.8±0.01 

0.21 

-0.75±0.14 

N/A 

N/A 

5.4 

-0.75 ±0.14 

-0.8 

4.56127 


omtff 

-0.77±0.1 

0.17 

-0.65±0.06 

N/A 

N/A 

7.1 

-0.8 ±0.0 

-1.0 

5.28883 


O.lbtff 

-0.38±0.04 

0.15 

-0.48±0.12 

N/A 

N/A 

30 

-0.75 ±0.03 

-0.6 

6.73403 


omtff 

-0.16±0.02 

0.04 

0.1±0.13 

176 

0.164 

206 

0.85 ±0.01 

-0.9 

5.46192 


omtff 

0.78±0.1 

0.03 

0.65±0.48 

350 

0.083 

655 

0.89 ±0.0 

-0.0 

4.45580 


0A5tff 

0.64±0.1 

0.05 

0.79±0.2 

2945 

0.010 

784 

0.89 ±0.0 

-0.7 

5.62422 


omtff 

0.99±0.09 

0.16 

0.64±0.18 

441 

0.066 

1175 

0.89 ±0.0 

-0.7 

12.7582 


omtff 

0.91±0.03 

0.24 

0.79±0.19 

5434 

0.005 

1939 

0.89 ±0.0 

-0.5 

8.81122 


trum of gravoturbulence using the t=0 snapshot of the 
column density (i.e. the timestep which has no gravity) 
and the maximum pixel value of the snapshot at some 
later timestep. We will call this map the “to + max(t)” 
column density map as it is the original t=0 map with 
one pixel added to it that represents the largest value of 
a snapshot at a later time. 

Figure [TT] plots the power spectrum similar to Figure 
M but overplots the power spectrum of the to+max(t) 
column density map. In general, the correspondence is 
fairly good for the both the later timesteps and and ear¬ 
lier time steps but has trouble matching the intermedi¬ 
ate timesteps where the turnover from turbulence domi¬ 
nated spectra to gravity dominated spectra takes place. 
At these time steps the power spectrum has a more com¬ 
plex behavior, i.e. can not be fit with a simple power law, 
which is evident by the lower value of kc in Equation [2TJ 

We fit the to+max(t) column density map using Equa¬ 
tion!^ and report the slopes (denoted as 13s) along with 
the values of the maximum of the column density map in 
Table H As is evident from visual inspection of Eigure 
El the largest discrepancies between the actual column 
density map and the to-l-max(t) column density map are 
at the intermediate time steps (t=0.25t//-0.35t//). Ad¬ 
ditionally at these time-steps the maximum value of the 
column density map (Table [2l column 7) increases dra¬ 
matically, signaling that gravitational collapse is begin¬ 
ning to dominate the low k values of the power spectra of 
the gas. We note that our results would be similar if we 
had used several delta function-like points or if we had 
placed them in different areas of the map. 

The fact that one can alter the turbulence power spec¬ 
trum to mimic the gravoturbulence power spectrum pro¬ 
vides researchers with an avenue of comparing the power 


spectrum of turbulence simulations with self-gravitating 
clouds following the prescription of maximum values 
given in Table [21 It also indicates that, on the scales 
of GMCs, the formation of cores introduce (5-function 
like intensity profiles to the power spectrum. Einally, if 
it is possible to mimic the effects of the gravoturbulence 
power spectrum then it should be feasible to remove the 
signatures of gravity as well. We discuss this in the next 
subsection. 

4.3. Restoring the Turbulence Power Spectrum in 
Self-gravitating Clouds 

In Section we demonstrated that it is possible to 
reproduce the self-gravitating power spectrum from the 
turbulent power spectrum, by adding a suitably chosen 
S function to the turbulent spectrum. This suggests that 
the removal of a (5 function from the self-gravitating spec¬ 
trum can recover the turb u lent s pectrum. It was first 
shown by iBeresnvak et al.l l)2005h that, for supersonic 
turbulence, the power the spectrum of p is made steeper 
by reducing rare density peaks, either by restriction or 
the logarithm. Here we try to use the relatively flat col¬ 
umn density spectrum to recover the steeper turbulent 
by a suitably chosen restriction. 

We define 


to be truncated column densities, and P<c{k), P>c{k) be 
the power spectra of each. We saw in Section 14.21 that 


S < Ec 

E > Ec 

E < Ec 
E > Ec 


( 22 ) 

(23) 
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one could construct a new field, 

El = Ensg(x) + S max (t)(5(x) (24) 

where Snsg is the column density before the action of self 
gravity, and Ei„ax(t) is the peak column density at some 
time t, then the power spectrum of Ei is quite similar to 
the power spectrum of the self gravitating cloud at time 
t. Since it can be shown that the power spectrum of i5(x) 
is a constant with respect to fc, we aim to select a cutoff 
column density, Ec, that has a power spectrum that is 
similarly constant over a given range. Thus for some self- 
gravitating cloud, if P>c{k) = const for some threshold 
Ec, then we show that the spectrum of the lower column 
density gas, P<c recovers the initial turbulent spectrum 
reasonably well. We show this in Figure [T^l In this 
figure, the left column is the high magnetic field simula¬ 
tion, and the right is the low magnetic field simulation. 
The top row shows P>c, which c is chosen to produce 
flat spectra; the bottom row shows P<C, the truncated 
spectra. The black line is Ps^sg, the initial turbulent 
state, and the colored lines are spectra for subsequent 
collapse states, as in earlier figures. To ensure P>c is flat 
we used a recursive bisection technique, varying Ec until 
the slope of P>c for fc > 10 was near zero. While there 
is significant evolution in the full column density power 
spectrum (e.g. Figure [9]), the evolution of the slope of 
the truncated spectrum is significantly restricted. While 
this technique is promising, it is is not yet able to dis¬ 
criminate between the initial turbulent conditions. The 
slope of the full spectrum, P{k), for the turbulent state in 
the high field run, for \ogk between 1.1 and 1.4, is —0.5, 
while the truncated slopes vary from —0.8 to —0.2. For 
the mid field, the full slope is —1.0 and the variation is 
from —1.0 to 0.1. For the low field, the full slope is —0.8, 
and the truncated slope varies from —0.9 to 0.0. The 
last column of Table [2] shows the fits. We hope that with 
some refinement this technique will allow for the further 
separation of collapsing and turbulent clouds. 

5. DISCUSSION 

Turbulence in GMCs is believed to be a part of a 
casca de that extends over 12 orders of magnitude in 
scale (|Armstrong et al.l 1199511 . In this paper we study 
the observable signatures of self-gravitating magnetized 
supersonic turbulence by applying the probability den¬ 
sity functions (PDFs) and the spatial density power 
spectrum to synthetic column density maps generated 
with the Enzo code. Unli ke other recent studies (e.g. 
iFederrath fc KlessenI (|2013lU that relate density diagnos¬ 
tics to star formation efficiencies/rates that rely on sink 
particle prescriptions, we stick entirely to observable col¬ 
umn density statistics and do not use any numerical pre¬ 
scriptions for star formation. In this sense we study the 
effects of gravitational collapse on driven turbulence and 
the evolution of the cloud as collapse proceeds. 

We find that there exists three characterizable stages 
of the evolution of the PDFs and density power spec¬ 
trum of the collapsing cloud which we term ” early,” ” in¬ 
termediate,” and ’’advanced.” The natural question that 
arises is how well do the predictions made in this pa¬ 
per match up with observations? The PDF of molecular 
clouds has been studied in a number of recent surveys 
and we have over-plotted mea surements of the powerl aw 
tail index of those published in iSchneider et al.l (|2014al lbll 


Figure [31 These clouds include NGC2603, which is a 
high-mass active star forming regions, the Auriga cloud, 
which is a low-mass star formation region, and the Orion 
B and Aquila clouds, w hich are moderat ely star forming 
( Schneider et al.l l201,3[l . Additionally, iSchneider et al.! 
(1201.311 published the Herschel data for the low-mass star 
forming cloud called Maddalena, which they found has 
a power law slope values of -3.65, however we do not in¬ 
clude this value in Figure|3]as it is steeper than the range 
plotted. 

How well do the ages and star formation history match 
with trends predicted by the simulations in Figure |3|? We 
discuss some of the literature on each of these clouds in 
order of evolutionary age predicted by the PDF in Figure 
[3] (youngest to oldest). 

1 . The Maddalena cloud has 41 young stars with disks 
and 33 protostars in the cen ter of the cloud with an 
age estimate of a few Myr (jMegeath et al.ll200^ . 

2 . Aquil a’s age estimates are expected to be a few 
Myr (jPrato et al.ll2008l l. 

3. Auriga-Galifornia age estimates are not well con¬ 
strained but it is suggested this cloud is not 
very evolved based on a high fraction of Glass I 
and Glass F YSOs (iBroekhoven-Fiene et al.ll2014f) . 
iHarvev et al.l ()201.3[ ) tab ulated 60 likely pr e-main- 
sequence objects while (jLada et al.l l201(lll report 
175-279 YSOs. 

4. Orion B has approximately 635 YSOs (iLada et al.l 
[2013). 

5. The age estirnates o f NGG3603 are 10 — 20Myrs 
(jBeccari et al.l I2010H making it the oldest cloud 
of the ones we compare with here. There are 
more than 10,000 stars with 7000 young stars 

forming a single power la w IMF) in NGG3603 
Haravama fc Martin3l2008il . 

In light of the literature on these clouds it would seem 
that the PDFs are able to trace the evolutionary state of 
the cloud. This bodes well for future observational stud¬ 
ies, which should combine the measures presented in this 
paper including the PDF variance, skewness, kurtosis, 
Tsallis and power spectrum to dissect the evolutionary 
state of clouds. 

Our results also revisit the now well established debate 
over the impact of magnetic fields on the star formation 
process. The idea that magnetic fields wholly dominate 
the formation of stars is accepted with less enthusiasm 
these days. Indeed, within the dynamical picture of the 
interstellar medium there are factors in addition to the 
magnetic field that impede the gravitational collapse of 
a molecular cloud and make star for mation less efficient, 
comp arable to the observed rates Isee lMcKee fc Ostriked 
120071) . For instance, turbulent magnetic fields u ndergo- 
ing fast reconnection (jLazarian fc Vish niaclll999il. whic h 
results in reconnection diffusion (jLazarian et al.l [20121) . 
diffuse out of clouds orders of magnitude faster than the 
typical ambipolar diffusion timescale that is usually in¬ 
voked in the traditional theory of molecular cloud evolu¬ 
tion. Our results here show that it is wrong to disregard 
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Figure 12. Top row: “upper” spectra, P>c(A:) which has been constructed to be flat. Bottom row: “truncated” spectra, P<c(fc). Most 
of the evolution of the cloud due to the collapsing gas has been removed in this bottom row. Left column shows the high field run, while 
the right column shows the low held run. 


the role of magnetic fields completely, as we see impor¬ 
tant dependencies in the evolution of the density statis¬ 
tics on the magnetization of the media. Furthermore 
we note that Table 2 shows values of our dense core re¬ 
gions being within the range predicted for gravitational 
free fall to become faster than the reco nnection diffu¬ 
sion ti mescale (i.e. around see lLazarian et aP 

(j2012l) . Future works should quantify the role of recon¬ 
nection diffusion in gravo-turbulent simulations such as 
the ones studied in this work. 

Currently the gravitational state of a molecular cloud 
is found by comparing kinetic and gravitational energy 
through the virial parameter, a = ba'^R/M, which re¬ 
quires knowledge of the clouds linewidth, cr, size, R, and 
mass, M, as well as the implicit assumption that a and 
M are spatially correlated. The techniques presented 
here provide a potential tool to probe the evolutionary 
state of a cloud using only the column density, which 
eliminates several sources of uncertainty in the estima¬ 
tion and implicitly includes information about magnetic 
field strength. Future numerical studies should focus on 
varying the sonic and Alfven Mach numbers in order to 
change both the power spectral slope and dynamical im¬ 
portance of the magnetic field in order to determine how 
turbulence speeds up or impedes star formation. The 
filtering techniques discussed in this paper can illumi¬ 
nate the dynamics of turbulence even if the cloud is star 
forming and in an advanced state of collapse. 

The current simulations aim to isolate the effects 
of magnetic fields and gravity on supersonic turbu¬ 


lence. However the use of periodic boundary condi¬ 
tions, solenoidal driving, and the lack of feedback are 
potential sources of discrepancy with real star forming 
clouds. Real clouds form from some sort of compres¬ 
sive flow, be it cloud coll isions, gravitation al instabil¬ 
ity, or thermal instability (iDobbs et al.ll201^ and refer¬ 
ences therein) and this has been shown to impact powe r 
spectra in such simulations (jFederrath fc Klessenll201,‘lh . 
However it is possible that in the inertial range the 
solenoidal-to-compre ssive ratio reaches a universal value 
(iKritsuk et al.l 1201011 . Further the periodic boundary 
condition imposes a conservation of total volume that 
may artificially enhance the amount of low density gas 
in the cloud. The impact of feedback is in general to 
inject kinetic energy, which st eepens the slope by sup¬ 
pressing small scale structure (|Sun et al.l 1200611 . Future 
simulations will incorporate these effects. 

Although the literature is abounding with papers re¬ 
garding the 3D density PDF and its relation to the star 
formation efficiency, the statistical properties of velocity 
and magnetic field are also of vital importance to describe 
and understand the turbulence cascade and the star for¬ 
mation process. Further, these statistical properties are 
essential to differentiate between the ’’log jam” of the¬ 
oretical models of star formation. Much work has gone 
into developing a toolbox of measurements to ascertain 
the physical conditions in GMCs and the ISM. The work 
presented here is a step towards extending this toolbox to 
self-gravitating turbulence. Techn iques such as th e Ve¬ 
locity Coordinate Spectrum fVCS. iLazarianl (|2009ll l can 
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provide the injection scale of turbulence and turbulence 
energy density and compare the observed spectrum of 
clouds to analytical predictions. The studies presented 
here can be improved upon in the future using tech¬ 
niques that go beyond the PDFs and the power spectrum; 
for instance, earlier research with non-self-gravitating 
turbulence has shown that the anisotropy of velocity 
fluctuations can be used to find the magnetization of 
the medium (iLa zarian et al.l 120011: lEsauivel &: LazarianI 
1200,1120101 : iBnrkhart et al.ll2014H . especially if used in 
conjunction tools to m easure the sonic Mach number (see 
IBnrkhart et al.ll201.1l and references therin). 

It is worth summarizing the impact of the three pa¬ 
rameters explored here (Mach number M.s: Alfven Mach 
number M.a, and collapse state) on each of the diagnos¬ 
tics available to us: 

1. Column Density PDF power law tails depend 
strongly on collapse state and Ma- 

2. Column Density variance from log-normal fit de¬ 
pends strongly on Ms, slightly on Ma, and little 
on collapse state. 

3. Column Density variance from direct measurement 
of column density is sensitive to the collapse state 
and Ms, and moderately correlated with Ma- 

4. Skewness and Kurtosis of f are not strongly cor¬ 
related with Ms', they are not strongly correlated 
with Ma, but are strongly correlated with collapse 
state. 

5. Tsallis fit parameters: q, w, A depend strongly on 
collapse state and Ma- 

6 . Density Power Speetrum is a very strong indicator 
of collapse state. 

7. Velocity Power Spectrum: While not stud¬ 

ied in this work and inclu ded here for com¬ 
plete ness, it was shown bv (I Collins et ahl l2012f) 
and l|Federrath fc Klessenll201^ that the velocity 
power spectrum is sensitive to A4s and Ma, but 
not the collapse state. We believe that the observed 
low correlation of the velocity power spectrum with 
the collapse state is the consequence of the rela¬ 
tively low energy injection to the turbulence cas¬ 
cade during the collapse as compared to the driving 
energy. This matter should be explored further in 
the future in simulations and in observat ions using 
the VCS/VCA methods (|Lazarianll2009ll . 

The differential sensitivity of the column density statis¬ 
tics on the key parameters of self-gravitating turbulence 
opens up prospects for determining the sonic Mach num¬ 
ber, Alfvenic Mach number and the collapse state by 
combining different statistics. 

6. CONCLUSIONS 

Turbulence, magnetic fields, and gravity are some of 
the key ingredients in the star formation process. Using 
synthetic column density maps generated with the Enzo 
AMR code, we investigated the robustness of the proba¬ 
bility density functions and the spatial power spectrum 
for understanding and separating the roles of magnetic 


fields, supersonic motions and gravitational collapse in 
the observable column density distribution. The PDFs 
and power spectrum reveal three stages of cloud evolu¬ 
tion as it progresses from diffuse turbulence dominated 
to collapse dominated. 

Regarding the PDFs, the statistical moments and the 
Tsallis fit to the observable column density distribution 
we find that: 

1. For early times, i.e. t < O.lSt//, the cloud has a 
lognormal distribution and will develop a high den¬ 
sity power law tail at intermediate times between 
t=Q.25tff — 0.35t//. The development of the tail, 
including its slope, also depends on the magnetiza¬ 
tion of the cloud. The PDF power law tail forms 
earlier in clouds with lower magnetization. The 
tails then become increasingly shallow as the col¬ 
lapse proceeds with time. 

2. The directly calculated variance of the column den¬ 
sity map (and natural logarithm of the column den¬ 
sity map), is a sensitive diagnostic for the cloud 
evolution. The variance increases monotonically 
with time as the cloud collapses and depends on 
the magnetic held strength. This increase in vari¬ 
ance with collapse causes the cloud to deviate from 
the sonic Mach number - variance relationship ex¬ 
pected from non-self-gravitating turbulence. 

3. The skewness and kurtosis of the natural logarithm 
of the column density map are insensitive to the 
sonic Mach number and trace the collapse of the 
cloud with time, i.e. they track the formation of 
the power law tail. 

4. The three Tsallis Ht parameters for the incremen¬ 
tal PDFs all strongly trace the evolution of the col¬ 
lapse with time: higher values correspond to more 
evolved cloud collapse. A strong magnetic field in a 
gravoturbulent cloud produces lower values of the 
amplitude and width of the incremental PDFs. 

We find the spatial power spectrum to be complimen¬ 
tary to the PDFs for studies of the evolutionary state of 
collapsing clouds. In particular we find that: 

1. The column density power spectrum of supersonic 
self-gravitating turbulence shows characteristics of 
a turbulence only power spectral slope at early 
stages of collapse (rs O.lSt//). At intermediate 
time-steps t=0.25t// — 0.35t//, the inertial range 
slope becomes increasingly shallow and the dissipa¬ 
tion range curves upwards. Eventually at advanced 
times the slope becomes positive and a positive 
sloped power law is seen down to the dissipation 
scales. The timescales of the changes in the slope 
depend on the magnetic held, with lower magnetic 
held facilitating earlier increases in the power spec¬ 
tral slope. 

2. We ht a three parameter function to the gravo- 
turbulence power spectrum: P{k) = Aik^'^e~^/^‘^, 
where Ax describes the amplitude, describes 
the power law behavior and the scale kc character¬ 
izes the intermediate stage turn over from turbu¬ 
lence dominated to self-gravity dominated which is 
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observed in the power spectrum at times between 
« 0.25tff — 0.35tff depending on the strength of 
the magnetic field. The exponential term only be¬ 
comes important (i.e. kc becomes on the order of 
100 ) when the first cores begin to form in our sim¬ 
ulations. The power law slope of the fit, /32 is com¬ 
parable to the traditional linear fit power law slope 
in the inertial range. 

3. We find that the effects of self-gravity on the power 
spectrum can be mimicked in a turbulence only 
simulation by including a single point in the col¬ 
umn density map that is equivalent to the maxi¬ 
mum valued point in a self-gravitating map at a 
given time-step. This addition of a point source 
is reminiscent of including delta-function like be¬ 
havior to the power spectrum at high k values. We 
provide values of the maximum points at each given 
time step and find that the slopes reasonably match 
the time-steps when kc is very large and the expo¬ 
nential term to the gravoturbulence power spectra 
relation is negligible. 

4. We find that the effects of self-gravity on the power 
spectrum can be removed and the turbulence power 
spectrum restored through spatial filtering of the 
high density material. 
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